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ABSTRACT 

Methods for simulating dependent sequences of contin- 
uous positive-valued random variables with exponential 
uniform. Gamma, and mixed exponential marginal distributions 
are given. In most cases the sequences are first-order, 
linear autoregressive, Markovian processes. A very broad 
two-parameter family of this type, GNEAR(l), with exponential 
marginals and both positive and negative correlation is 
defined and its transformation to a similar multiplicative 
process with uniform marginals is given. It is shown that 
for a subclass of this two-parameter family extension to mixed 
exponential marginals is possible, giving a model of broad 
applicability for analyzing data and modelling stochastic sys- 
tems, although negative correlation is more difficult to obtain 
than in the exponential case. Finally, two schemes for auto- 
regressive sequences with Gamma distributed marginals are out- 
lined. Efficient simulation of some of these schemes is discussed. 
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1 . 



INTRODUCTION 



In a recent series of papers 11,2,3,4,5,6,7,8,9] some 
simple models have been derived for stationary dependent 
sequences of positive, continuous random variables with given 
first-order marginal distributions. In general the dependency 
structure, as measured by second-order joint moments (serial 
correlations) mimics that of the usual linear mixed auto- 
regressive-moving average (ARMA) models which have been used 
for so long in time-series analysis. In the ARMA models, 
which are defined quite generally, there is in usage an 
implicit assumption of marginal normality of the random vari- 
ables. This is clearly not the case if the random variables 
are positive, say the times between events in a series of 
events (Cox and Lewis [10]) or the successive response times 
at a computer terminal. Thus the new models are derived to 
accommodate situations in which the dependent random variables 
have, for instance, exponential, Gamma, uniform and mixed 
exponential marginal distributions. The exponential case is 
the most highly developed, with the nomenclature (Lawrance and 
Lewis [4]) EARMA (p,q) (exponential process with mixed moving 
average-autoregressive structures of orders p and q respec- 
tively) and NEARMA ( p , q) (new EARMA(p,q)). A generalization 
to extend the range of attainable autocorrelations to nega- 
tive values has been defined with the nomenclature GNEARMA(p,q) . 

The development of the probabilistic properties of 
these processes is given in the referenced papers, applications 
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to queueing models and computer systems modelling by Lewis 
and Shedler [11] and Jacobs [12,13], while development of 
estimation and testing procedures has just begun. 

The object of the present paper is to define and 
discuss the simulation of the processes on digital computers, 
though for the sake of brevity only the first-order Markovian, 
autoregressive case is considered. The simplicity of structure 
of these models — in general they are linear additive mixtures 
of random variables — makes them ideal for this purpose. How- 
ever, stationarity conditions are sometimes difficult to derive 
analytically and in some cases it is not simple to generate 
the innovation random variables in the processes. A striking 
example of this is the case of the Gamma first-order autore- 
gressive process, given in Section 4A, for which an efficient 
means of simulation was reported by Lawrance [7] for some para- 
metric values. This procedure carries over into another Gamma 
process, the Gamma-Beta process, which will be discussed in 
Section 4B. 

In Section 3 it is shown that a simple transformation 
of the exponential sequences gives a direct multiplicative 
method for generating dependent processes with uniform marginals. 
These could be the basis in simulations for many other types 
of dependent sequences. 

Finally the NMEAR(.l) process is detailed in Section 4C; 
this generates a first-order Markovian process with mixed expo- 
nential marginals. It is useful for simulating situations in 
which the observed random variables are correlated and over- 
dispersed relative to an exponentially distributed random variable. 
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2 . 



EXPONENTIAL AUTOREGRESSIVE MARKOVIAN SEQUENCES 



We give here three methods of generating first-order 
autoregressive, Markovian sequences with exponential marginal 
distributions. The first two are defective in terms of their 
sample path properties (.the first more so than the second) 
while the third, NEAR(l) and its generalization GNEAR(l) , is 
satisfactory in this respect and is a very rich model. The 
defect of the first two sequences is also highlighted by the 
simulation procedures used; they can be generated from one 
sequence of exponential variables. 

The word "autoregression" in the context of a stochas- 
tic sequence {X n } is often used rather vaguely. In the 
first place linear, additive autoregression is usually implied. 
In the second place first-order autoregression can mean that 
in the defining equation for X r the previous value enters 
explicitly. Thirdly, it can mean that the conditional expec- 
tation of X , given X . = x ,, is an additive linear 
function of x n _]/ 



E ( X |X , = x .) = 
n 1 n-1 n-1 



a + bx 



n-1* 



( 2 . 1 ) 



The processes discussed in this paper are autoregressive in 

the latter two senses and, except in the case of uniform 

marginals, are autoregressive in a linear additive way. They 

are also Markovian; the Markovian property (first-order) means 

that the probability structure of X n , X n+ ^,... > given 

X^ . = x„ is independent of X ~ , X„ , . . . . 

n-1 n c n-2 n-o 
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2A. The Exponential DAR(l) Process 

A very simple exponential autoregressive Markovian 
sequence is generated by the equation (Jacobs and Lewis 
[14,15]) 



n 



= V X i 
n n-1 



+ (1 - 



V )E 
n n 



( 2 . 2 ) 



where the V 's, n = 1,2,... are i.i.d with 
n 

P{V n =l} = 1 - P{V n =0} = p and E n , n = 1,2,... are, as 
throughout the paper, independent exponential random variables 
with parameter A and independent of the ' s ; that is 



P{E < x} = 
n — 



1 - e 



-Ax 



x > 0 , 



A > 0 



(2.3) 



= 0 , x < 0 . 



For this process the serial correlations = corr (X fi , X^ +k ) 

are 

P k = P k (2.4) 

and 

E(X nl X n-l = x n -l } = p l X n-l + (1 " P l )/A ‘ (2 ‘ 5) 

This process, which was introduced to model discrete valued 

variables, is not well suited to modelling continuous data 

because runs of X 's with the same value can occur quite 

n 

frequently in the sample paths of the process. This happens 
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when x n _i is picked successively in (2.2), rather than the 
innovation E n> Moreover the lengths of the runs of similar 
values are geometrically distributed. 

2B. The Exponential EAR(l) Process 

Another model is derived from the usual linear model 



X 



n 



pX 



n-1 



+ 



e 

n 



( 2 . 6 ) 



in which the i.i.d. innovation process {e } is chosen so 

n 

that the x n ' s are marginally exponential (A) . Gaver and 
Lewis [1] show that for this to be true, one must have 
0 < p < 1 and 



e 



n 




w.p. 1-p, 



w.p. p , 



(2.7) 



where {E n )r as previously, are i.i.d. exponential (A) . 

Again p^ = p and E ( x n l x n _i = x n _i) = p l x n -l + ( 1_p i)^' 
as at (2.4) and (2.5) for the exponential DAR(l) model. The 
difference is in the sample paths; the EAR ( 1 ) process simula- 
tions show runs of geometrically decreasing X r ' s , but no runs 
of constant value. The geometrically distributed runs occur 
when only pX n _p is picked in (2.6). 

The Markov property of the two sequences implies that 
if Xq is chosen to be Eg, an exponential ( A) random variable 
independent of E^, E 2 , ... , then X^, X 2 , • • • forms a 
stationary sequence. 
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Naive inspection of the defining equations (2.2), (2.6) 



and (2.7) suggest that to generate a stationary sequence of 
length N, X^..., X^ , (N+l) i.i.d. exponential deviates and 

N uniform variates (for the selection process) are needed. 
However, the sequences can be generated from only one exponen- 
tial sequence; this is possibly related to the degeneracy in 
the processes. This method uses the memoryless property of 
exponential (A) variables, namely that if E n is given to be 
greater than a constant y, then E n - y is again exponential ( A) 

Thus the algorithm for the EAR(l) process is to ini- 
tialize by setting X Q = E Q ; subsequently set X n = pX n _ 1 if 

E < x = -£n(l-p)/A; otherwise set X = pX . + (E - x ) . 
n — p n n— i n p 

This uses the fact that, from (2.3) , P{E £ x } = p. 

n p 

Even greater efficiency can be obtained, though this 
must be qualified by considerations as to whether (i) the 
X n 's are to be generated one at a time or in an array; (ii) 
a subroutine is available to generate exponential random 
variables faster than can be done by taking logarithms of 
uniform deviates, and (iii) the speed of division in the 
computer is short compared to the time needed for generation 
of uniform deviates. 

The more efficient scheme recycles uniform variables, 
i.e. if U is given to be between constants a and b, 
where 0 < a < b _< 1, then (U-a)/(b-a) is a uniform random 
variable. (Note that its value is not given, only that it is 
in (a,b) ) . The expected number of uniform deviates required 
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to generate an EAR ( 1 ) process of length N with this 
algorithm is 1 + (l-p)N, which is less than the number N 
required to generate an i.i.d. exponential (1) sequence . Also 
the expected number of logarithms is (l-p)N, while N 
comparisons are always needed. 



2C. The Exponential NEAR(l) Process 

A broader two-parameter exponential sequence which 
is a first-order autoregressive, Markovian process and an 
additive linear mixture of random variables is given by 
Lawrance [7] and developed by Lawrance and Lewis [5]. Called 
NEAR ( 1) , the sequence is defined as 



X = 
n 



ex 



n-l 



e + s 
n 



w.p. 



w.p. 



a 



1-a 



n = 1,2 



/ • • • 



( 2 . 8 ) 



where 0 £ a _< 1 and 0 _< 3 1 but a = 3 =f 1. Also the 

selection process is done independently for each n. It can 
be shown that for the X n to be marginally exponential ( X) 
the innovation variable e n must be generated from an E n by 
the exponential mixture 



e n = 



n 



w.p. 6 = 



1-3 



1 - (1 - a) 3 



n = 1,2, 



(2.9) 



a3 



(l-a)3E n w.p. 1-6 = YZ J i _ 



a) 3 



8 



providing a and 3 are not both equal to one. When a = 0 
or 3=0 the {X n > are i.i*d* exponential variables, 
whereas when a = 1 the EAR(,1) model given at (2.6) and (2.7) 
is obtained. In fact choosing a as a function of 3 in a 
suitable way, e.g. 3 = a, gives an exponential model with a 
full positive range of serial correlation of order one, since 
it is easily shown that 

P k = (a3) k . (2.10) 

Again 

E ( X | X -| = x -, ) = a3x , + U-a3)/A 
n 1 n-1 n— 1 n- 1 

= Pi x n _i + (l-P 1 l/A 2.11) 

and Xq = E q gives a stationary sequence. Thus the correla- 
tions and regressions are the same as for the first two models. 
However the NEAR (.1) process allows one to model a broader class 
of exponential sequences, as measured either by sample path 
behavior or higher-order joint moments; see Lawrance and Lewis 
[5] for details. 

A particularly simple case occurs when 3 = 1; this 
model, called TEAR(l) , is very tractable analytically and, as 
will be shown in Section 4C, extends easily to the case of 
mixed exponential distributions for the X^. 

Note that in the NEAR(l) process the innovation e n 
is always present unless a = 1 and it is therefore not 
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possible to simulate the stationary process with less than 
N+l uniform variates. A detailed algorithm is given in 
Section 2E for a more general case. Since for a stationary 
array of N X R ' s , exactly N+l uniform deviates are required 
because of the ability to recycle the uniforms and transform 
them into exponentials, it could be advantageous to generate 
these uniform deviates in an array which would be replaced 
one at a time by the X 's. Care must be taken with the re- 
cycling of the uniform variates U if y = 1-a is close to 
one or zero. In that case it is probably better for computa- 
tional reasons to use 2 (N+l) uniform variates. Note that 
a = 0 gives the EAR(l) process. When 3=1 a simpler 
algorithm can be used since E n is no longer a mixture of 
two exponentials; this important special case is called the 
TEAR (1) exponential process. 

When 3 = 2 ~ another one-parameter sub-class of the 
NEAR(l) process is obtained with strikingly regular sample 
path properties. For this process P{X n > x n _i^ = 1/2; this 
is in striking contrast to the EAR(l) process whose sample 
paths show runs-down and the TEAR(l) process whose sample 
paths show runs-up. This is illustrated in Figure 1; all 
sequences have p = 0.75 and are transformations of the 
same E n sequence. The parameter space of the NEAR(l) process 
is illustrated in Figure 2. 
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ALP HA = .866 BETA= .866 ALPHA= .857 BETA= .875 



Figure 1. Sample paths of NEAR(l) processes, all with 
a£ = p 1 = 0.75, for different values of a and 6. Figure A 
is the 1 EAR(l) process (a = 1, 3 = 0.75), Figure B. is the 
TEAR ( 1 ) process (a = 0.75, 3=1)/ Figure C is the PREAR(l) 
process a = = 0.857, and Figure D is the REAR ( 1 ) process 

a = 3 = (0.75) x / 2 . All the sample paths are transformations 
of the same i.i.d. exponential sequence E^ , n = 0,1,..., 100. 
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0 



0 = 1 ; TEAR ( 1 ) 




a 



Figure 2. Parameter space for the NEAR(l) exponential auto- 
regressive process. The cases a = 0 and/or 6=0 give 
i.i.d. exponential sequences while B = l/(2-a) gives the 
partially reversible PREAR(l) process for which P(X n < X n _^) = 1/2. 
For a = 1 we have the original exponential process EAR(l) 
which tends to have runs-down. The TEAR(l) case, 3 = 1 
gives very simple analytics but exhibits runs up. Also shown 
is a locus of constant , in this case = aB = .5 . 
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2D. The Generalized Exponential Process GNEAR(l) with Negative 
Correlation 

The exponential processes defined in Sections 2A, 2B , 
and 2C do not exhibit negative correlation or alternation of 
correlations. Such behaviour is found in, say, a normal linear 
first-order process for which = cor (X n , ^n+j^ = and 

-1 < p < 1 so that, for instance, p^ can be negative. A 
scheme for broadening the correlation structure of the EAR(l) 
process is given in Gaver and Lewis [1], However in the ex- 
ponential case a much simpler alternative method is available. 

Assume for simplicity that A = 1. Now X , is 

n— 1 

a unit exponential variable and U n = F(X n _^) = 1 - exp(-X n _^) 

is a uniform (0,1) variable, as is 1 - U = exp(-X ,). Then 

n n- 1 

X^_ 1 = F _1 ( 1 - U n ) = F~ 1 ( 1 - F(X n _ 1 )) = - Zn ( 1 - exp(-X n _ 1 )) is 
a unit exponential variable; in fact it is the antithetic 
of x n _i which gives the maximum negative correlation 
attainable in a bivariate exponential distribution: 



r = corr(X . , X ,)=1 -tt/6 = - .6449 
n-1 n-1 ' 



(2 



Now the process 



• -X 

X =fe +BX n =e - 3 £n(l - e n ) 
n n n-1 n 



n 



w.p. a 
w.p. 1-a 



(2 



or 

X n 



e + U X , 
n n n-1 



(2 



. 12 ) 

.13) 

.14) 
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in which the U n 's are with P{U n =l} = 1 - P{U n = 0} = a 
gives a process with autocorrelations which alternate in sigh. 

In particular = r(a3). To combine this with the positive 
correlation case in a continuous way we introduce a new para- 
meter p e[0,l] and i.i.d. indicator variables I n , indepen- 
dent of U n , from which P{l n = 1} = 1 - P{l n = 0 } = P. Then 
the GNEAR(l) model is defined as 



and gives a complete range of first-order serial correlations 



Higher lag correlations are more complicated and will be given 

elsewhere. Note, however, that 1-p = l/(l-r) gives a case 

in which X and X . are a bivariate exponential pair which 
n n— 1 

are dependent but have zero correlation. 

Figure 3 give four sample paths for the case p = 0 
for values of a and 3 corresponding to those in Figure 1, 
which is the case p = 1. 

2e. Algorithms for the GNEAR(l) Process 

We give here two algorithms for the GNEAR(l) process, 
one based on generation of uniform deviates one at a time, 
the other based on the availability of subroutines to generate 
arrays of uniform and exponential variables. The sequences 
generated are unit exponentials (A = 1). The special case 
a = 1 (EAR ( 1) ) and 3=1 (TEAR ( 1) ) are handled separately as 
they will cause divides by zero in the main algorithm. 



X 



n 




(2.15) 



P 1 = a3 [ p + (1-p) r] . 



(2.16) 
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GNEAR(l) PROCESS 



P= 



.000 



A. 




ALPHA= 1.000 BETA= .750 




B . 




ALPHA = .750 BETA= 1.000 




Figure 3. Sample paths for the GNEAR(l) process with P = 0, 
corresponding to those in Figure 1 for the NEAR (1) process 
(p = 1). Here p x = - 0.6449 aB. 
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ALGORITHM GNEAR1A 



This algorithm generates a sample of size N from the GNEAR1 
process. It is based on the generation of uniform random numbers one at 
the time, with recycling of these uniform random numbers for further use. 

It is assumed that a subroutine (UNIFORM) exists that generates raw 
uniforms. Input values are N, ALPHA, BETA and P the last three taking 
values in the closed interval [0,1]. However ALPHA = BETA = 1 is not allowed. 



INPUT N, ALPHA, BETA, P 
CALL UNIFORM (U) 

X(0) <- -log e (U) 

C «- (1- ALPHA) *BETA 
IF ALPHA = 0 THEN 

ELSE 

END IF 

DO I «- 1 to N 

CALL UNIFORM (U) 

IF U < ALPHA THEN 



ELSE 



END IF 

IF U _< D THEN 

ELSE 

END IF 

X(I) «- Z + Y 



/* Generate uniform U*/ 

/* convert to exponential*/ 



D 1 

D (1-BETA) /(1-C) 



CALL UNIFORM (V) 

IF V < P THEN Y BETA* X(I-l) 
ELSE Y - BETA* log g 



1-exp (-X(I-l) 



END IF 

U •*- U/ALPHA /* Recycle U */ 

Y 0 

U «- (U-ALPHA) / (1- ALPHA) /* Recycle U */ 



Z «- -log e (U/D) 

Z •*- - C * log (U-D)/(l-D)) 



END DO 
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ALGORITHM GNEAR1B 



This algorithm generates a sample of size N from the GNEAR1 
process. It is based on the generation of arrays of uniform and expo- 
nentially distributed random numbers. It is assumed that subroutines 
UNIFORM and EXPON exist that generate arrays of uniform random numbers 
and exponential random numbers respectively. For each GNEAR1 number 
generated two raw uniform random numbers are used. The first is recycled 
in those cases where a third uniform is needed to make the selection with 
probability P. Input values are N, ALPHA, BETA and P the last three 
taking values in the closed interval [0,1], except the case ALPHA = BETA = 1 



INPUT N, ALPHA, BETA, P 
C (1 - ALPHA) *BETA 
IF ALPHA = 0 THEN 

ELSE 

CALL UNIFORM (U, 2*N) 

CALL EXPON (E, N+l) 

X(0) E (N+l) 

DO I -e 1 to N 

IF U(I) £ ALPHA THEN 



ELSE 

END IF 

IF U(I+I) £ D THEN 

ELSE 

END IF 

X(I) Z + Y 
END DO 



D + 1 

D (1-BETA) /(1-C) 

/*Generate array of 2N uniforms*/ 
/*Generate array of N+l exponentials*/ 



V = U(I) /ALPHA 
IF V £ P THEN Y + 
ELSE Y + 



/* Recyle U */ 



BETA* X(I-l) 



- BETA* log g 



1-exp (-X(I-l) ) 



END IF 



Y 0 



Z E (I) 

Z C*E (I) 
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3. UNIFORM MARKOVIAN SEQUENCES, NUAR(l) 

It is convenient to have dependent sequences of random 
variables with marginal distributions other than exponential. 

Before discussing other solutions to Equation (2.8) we show 
that a simple transformation of the NEAR(l) process gives a 
two-parameter family of Markovian random variables with 
uniform marginal distributions. It is well-known that an 
exponential transformation of a unit exponential random vari- 
able gives a uniformly distributed random variable. Thus we 
have from (2.8) and (2.9) the multiplicative model for a 
uniform Markovian sequence (X n ), n = 1,2,... , called 

NUAR(l) , 

( ft 

(3.1) 

(3.2) 
(3. 3) 

for U , n = 1,2,..., i.i.d. uniformly distributed, providing 
that a and 8 are not both equal to one. Again if Xq 
is uniformly distributed and independent of U^ ^ , ... the 
sequence is stationary. 



where 



X = 
n 



e a i 

n n-1 



n 



e = 
n 



U 



n 



U. 



( 1-a) 8 



n 



w. p. a 

w.p. (1-a) 
n _ 1,2,... , 



„ ^ v _ 1-8 
w,p ‘ Y “ l-(l-a) 8 



w.p. 1-Y = 



a8 

l-(l-a) 6 



n = 1,2,. 
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The sequence is clearly quite simply extended to give 

negative correlation, as in the GNEAR(l) process; in fact 

X n _^ in (3.1) is just replaced by (1-X n _^). Algorithms are 

easily obtained by adaptation of those given in Section 2E 

for the exponential case. It remains to find the correlation 

structure and the regression of X on X , . 

n n— 1 

To do the former, let X* be a NEAR(l) sequence with 
A = 1, so that the sequence X^ at (3.1) is given by 

r * 

X = exp{-X*}. Now the joint Laplace-S tielt jes transform of X 
n n n 

and X n _ k , <f> x * , * (s,t) = E{exp[- sX* - tX* ]}, is given 

n n-k 

by Lawrance and Lewis [5], Setting s = t = 1 in 

4> x * y* <s,t) gives, 
n' n-k 

^X* X* (Ifl) = E{exp(— X*) exp(— X*_^ ) 

= E(X n X n-k» • (3 - 4) 



Then using the fact that for a uniform random variable 

E(X) = 1/2 and var(X) = 1/12 we have, after simplification. 



p k = corr(X n ,X n _ k ) 



3 

2 + 6 k 



k 



1 

i = l 



l 1 + 



aB 

( 1-a) f3\ 



k = 1,2,... 



(3.5) 



Note that this is not simply a geometrically decaying corre- 
lation sequence, as for the NEAR(l) process. However, for 
the important special case when (3 = 1 we get 



19 



k 




a 



k = 1,2 



(3.6) 



2-a 



/ • • • 



and thus the serial correlations are the kth power of 

p^, which takes on any value between 0 and 1. Thus we have 
a particularly simple uniform Markovian sequence, although 
the sample paths will tend to have runs-up. 

A similar analysis given in Lawrance and Lewis {5] 
shows that 



so that the regression is not linear for this Markov process 
with uniform marginals, unless 3=1. 

This uniform sequence could form the basis, via a 
probability integral transform, of many other sequences with 



choice for a one parameter model since this case gives a sample 
path which is partially time-reversible (see Lawrance and 
Lewis [5]), with a balance of runs-up and runs-down. However, 
marginal transformations do not preserve correlation structure, 
as shown at (3.5) , and it is therefore useful to see whether 
sequences with marginals other than exponential can be generated 
from (2.8); this requires finding, if possible, a suitable 
choice of innovation sequence e n « The result will be a simple 
process with autoregressive Markovian structure and the desired 
marginal distribution. 




{1 - a + au^} 



(3.7) 



given marginals. The parametrization 8 = is a good 
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4. MARKOVIAN SEQUENCES WITH SOME OTHER MARGINALS 

Although an exponential distribution is a common assump- 
tion for positive random variables met with in problems in 
operations research, it is too narrow an assumption to encom- 
pass many real situations. Therefore parametric distribution 
models are invoked which include the exponential as a special 
case and which allow for the modelling of data which has greater 
or lesser dispersion than exponentially distributed data. Two 
commonly used models are 

(i) the Gamma (k, A) distribution whose probability density 
function is 



f(x) = 



, , , » k— 1 — Ax 
A ( Ax) e 

r(k) 



k>0; A>0;x>0, 



(4.1) 



where T(k) is the complete gamma function, and 
(ii) the (convex) mixture of exponential random variables 



-A,x ~^2 X 

f(x) = TT 1 A 1 e + (l-TT 1 )e , 0 < A 1 < A 2 ; 



X > 0, 0 _< TT^ £ 1 . 



(4.2) 



The Gamma distribution has dispersion, measured by the coef- 
ficient of variation C(X) = c(X/E(X)), which is greater than 
the exponential value of 1 if k < 0 and less than 1 if k > 1. 
The mixed exponential always has C(X) >_ 1 , the equality 
occurring when the special case of an exponential random 
variable with parameters A^ or A 2 holds. 
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4 A. The Gamma GAR(l) process 



The solution of the standard first-order autoregres- 
sive equation (2.6) with stationary gamma marginals defines 
the GAR(l) process. Using Laplace-S tielt jes transforms with 
(2.6) shows that for X n to be Gamma (k,A) , we must have 



4> e (s) = E (e~ S£ ) 




(1-P) 




(4.3) 



For k integer this has an explicit inverse. For example, 

2 

for k = 2 the innovation e is zero with probability p , 

is exponential ( A) with probability 2p(l-p) and is Gamma(2,A) 

2 

with probability (1-p) . It is easy to show in general that e 

is zero with probability p , so that the "zero defect" is not 
serious for large k. A method of simulating a random variable 
whose Laplace-S tielt jes transform is equation (4.3) was derived 
by Lawrance [7], using the fact that this sequence arises in 
a particular type of shot noise process. From this we have the 



Gamma Innovation Theorem 

Let N be a Poisson random variable with parameter 
8 = -k &n(p). Let U^, be uniformly distributed 

over (0,1) and independent. Let , . . . , Y N be exponential ( A) 
and independent. Then z can be simulated using 



e 






N 



l 

m=l 



U 



Y p 
nr 



m 



if N > 0 , 



= 0 



if N = 0 . 



(4.4) 
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A proof is not given here. Note that e is zero with prcb- 
ability exp{-k £n(p) } = p . Also the Poisson number N 
of uniform and exponential random variables which must be 
generated for each e has expected value 0 = -k £,n(p) . 

This will be prohibitively large, and the simulation will be 
very inefficient, if k is large and/or p is close to zero. 
Neither of these cases is serious, however. If k is large, 
say greater than 50, the sequence is almost normal and the 
usual normally distributed, AR(1) linear process can be used. 

If p is as small as 0.001 then E(N) is only k x (6.9078) 
which is still reasonable. However, for p this small the 
sequence is approximately i.i.d. Gamma and acceptance-rejection 
techniques for simulating Gamma variables are known. 

It is quite simple to write algorithms for the GAR ( 1 ) 
case analogous to those in Section 2E. It would pay to have a 
built-in routine for generating the Poisson variable which will 
bypass further calculations if N = 0. In other words routines 
for generating Poisson variates which start by searching at the 
median of a table of cumulative Poisson probabilities will be 
inefficient. 

Unfortunately the NEAR(l) process does not appear to 
extend to the Gamma case; it can be shown explicitly that 
there is no innovation e n in equation (2.8) which will make 
X have a Gamma distribution with k = 2 if a 4 1* 

There is, however, another model, the Gamma-Beta model, 
inspired by an example in Verwaart [16] and discovered inde- 
pendently by Fishman [17], which is quite broad and which can 
be simulated using the Gamma Innovation Theorem. 
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4B . The Gamma-Beta Model, GBAR(l) 



This model is a linear autoregressive process with 
random coefficients which includes the GAR(l) process but is 
of limited practical use in data analysis because its likeli- 
hood function is analytically untractable. Nevertheless it 
could be useful in simulations, particularly if the zero-defect 
in the GAR(l) process is unacceptable. 

Thus, we define the stationary sequence as 

X n = 6B n X n _ 1 +e n 0 < 3 < 1 , n = 0, +1, + 2,..., (4.5) 

where the B^s are i.i.d. and independent of x n _]/ and 

B , for X to have a Gamma (k,A) distribution, has a 
n n ' 

Beta(k-b,b) distribution with mean E (B n ) = (k-b)/k. The density is 

f B ^ “ r (k-b) r (b) X (1-x) . 0 < x £ 1 . (4.6) 

n 

0 < b _< k . 

The distribution of the i.i.d. e's to make X Gamma (k,A) 

is still to be determined. To do this and to see the rationale 

behind the model, it is simplest to obtain the distribution of 

BX n . 
n n-1 

Now B n can be generated as Z^/(Z^+Z 2 ) where 
Z^ and Z 2 are independent and Gamma (k-b, A) and Gamma (b,A) 
respectively. Moreover, B n is independent of (Z^+Z 2 ) , which 
could be used to generate X n _^. Then, 

B n X n-l = { z 1 /( z 1 +z 2 )} {Z i+ Z 2 } = Z ± is Gamma (k-b, A) . Using 
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this fact, which can also be shown analytically, and the 
defining equation (4.3) we have 



E 




♦x <s) 



n 



*B X^f't’e (s) 

n n-1 n 



(4.7) 



so that on using the fact that for a Gamma (k , A) variable the 
Laplace-S teilt jes transform is (A/A+s) , we have 



♦e <s> 



X 



n 



(s)/<f> B (3s) 

n n-1 




(4.8) 

(4.9) 



Thus e n is generated as the sum of a Gamma (b,A) variate 
and, from (4.3) , a Gamma innovation variable which can be 
generated from (4.4). 

For this model 



p . = corr (X X . . . 
3 n n+j) 






k-b 

k 



) 



j 



j — 0,1,2, 



(4.10) 



and 



E (X 

n 




1 



x 



n- 



1 



) 



P,x . + 
1 n-1 



(1— P 1 )E (X n ) . 



Also when b = 0 we have the GAR(l) model. If 3=1 the 
model requires only a Gamma (b,A) variate for £ n and a 
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in the simulation. The case 



Beta(k-b,b) variate for B n 
k = 1 gives an exponential process which is not the same as 
the NEAR(l) process. 

It should be noted that while the Gamma Innovation 
Theorem makes the Gamma-Beta model tractable from a simulation 
viewpoint, it is still difficult to unite down a likelihood 
function or to get negative correlation. Thus further 
developments are needed from the Gamma case. An algorithm is 
given for the Gamma Beta Model on the next page. 

4C. Mixed Exponential Markovian Processes MEAR(l) and TMEAR(l) 

In addition to Gamma processes, first-order autore- 
gressive Markovian processes with mixed exponential marginal 
distributions can be obtained from equations (2.8) and (2.9) 
in two special cases, and these sequences should be widely 
useful in modelling stochastic systems. 

(i) The case a = 1; MEAR(l) . 

In Gaver and Lewis [1] it is shown that the solution 
to the Laplace transform of e n for the linear model (2.6) 
is a constant p plus a (generally) non-convex mixture of 
three exponential functions. This can be shown to be a proper 
density function if p _< ^/^2' but it can also be shown that 
it is not a density function for all p less than one and 
greater than or equal to zero. More particularly, Lawrance [6] 
showed that unless A^ is much smaller than A 2 (and thus 
the X n are very over-dispersed relative to an exponential 
random variable) a solution exists for e n for all p. Thus 
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ALGORITHM GBAR1 (GAMMA-BETA MODEL) 



This algorithm generates a sample of size N from the GBAR1 
process using array generators POISSON, UNIFORM, EXPONENTIAL and GAMMA 
to produce Poisson, Real Uniform (0,1), Exponential and Gamma distributed 
random numbers. 

The following restrictions exist on the input parameters: 

0 _< B _< K 

0 < K 

0 < BETA < 1 



INPUT N, BETA, B, K 
TH -f- -(K-B) * log e (BETA) 
CALL POISSON (PSN,N,TH) 

CALL GAMMA (Z1,N, K-B) 

CALL GAMMA (Z2,N,B) 

CALL GAMMA (X(0) , 1,K) 

DO I 1 to N 

CALL UNIFORM(U,PSN (I) 
CALL EXPON(E,PSN(I)) 



/* 

/*Generate N 
/*Generate N 
/*Generate N 
/* Initialize 



/* 



/* 



Initialization */ 

poisson deviates with parameter TH */ 
Gammas with parameter K-B */ 

Gammas with parameter B */ 

X(0) */ 

Generate PSN(I) real (0,1) uniforms */ 
Generate PSN(I) Unit Exponentials */ 



Y -i- 0 

DO J 1 TO PSN(I) 

Y «- Y + E(J) * BETA ** U(J) 

END DO 

BN -s- Z1(I)/(Z1(I) + Z2 (I) ) /* Compute BETA Deviate */ 

X(I) BETA * BN * X(I-l) + Y 
END DO 
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we have a useful process, although again the zero-defect of 

order p is a problem. However, one case which cannot be 

simulated this way occurs when 1/^2 = 0 , i.e. X is zero 

with probability 1 -tt^, and exponential (A^) otherwise. This 

kind of situation occurs in practice as, e.g., the waiting time 

for an item in an inventory system. Fortunately it can be 

handled in the next case. 

(ii) The case 8 = 1; TMEAR(l). 

When 8=1 in equation (2.8) , a mixed exponential 

process TMEAR(l) is obtained which is extremely simple to 

simulate since the innovation e is just the mixture of two 

n 

exponentials for all 0 £ p < 1. Moreover, the process has 
no zero-defect. As discussed above, the sample paths will 
tend to "run up," but this is no great problem unless p is 
fairly large. Thus we have the following Theorem (Lawrance and 
Lewis [18]) which we state without proof: 

TMEAR(l) Theorem 

Let the first-order autoregressive, Markovian sequence 

{X } be defined by 
n J 

X = £ + V X , n = 1,2,3,... 

n n n n-1 ' ' ' 

where for the i.i.d. sequence V , n = 1,2,3,..., 

P{V =1} = 1 - P{V =0 } = a for 0 < a < 1. Then the sequence 
n n — 

(X n ) is stationary and has a (convex) mixed exponential 
marginal distribution with probability density function 
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f x (x) 



-A^x 

7T i A i e + 




x > 0 



(4.11) 



where 0 < A 1 < X 2 ; 0 < 7^ < 1 ; ^ + * = 1 u 1 and „ 1 

A 1 2 A 2 

if e n is i.i.d. and has a (convex) mixed exponential 
distribution given by 



-Y]* -Y 9 x 

f £ (x) = n 1 Y 1 e + n 2 Y 2 e x > 0 

with J 1 > y 2 > 0; r^, n 2 > 0 , ^ + n 2 = 1 



where 



p = E(X) = tt 2 P 2 + ^i^i ' 



b = M 1 + y 2 “ aM ' 

B = y l + y 2 “ M ; 
a = ( 1-a) y 1 vi 2 ; 

y 1 ,y 2 = ^k+^tb^a] > ; 

Y 0 = 7T 2^ j 1 + 7T 1^ j 2 '* 

"l = <y 1 -y 0 )/(y 1 -y 2 > 



"2 = <Y 2 -Y 0 )/(Y 2 -Y 1 ) 



(4.12) 

(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 

(4.18) 

(4.19) 



and Xq is independent of ... and has probability 

density function (4.11). 

Note that the special cases where tt 1 = 0 or tt 2 = 1 
give NEAR(l) exponential processes with parameters A 2 and 
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X 1 respectively. Thus they should be handled by Algorithm 2 
since they will cause computational problems. The case 
X = X also gives a NEAR ( 1) process and is excluded for 
similar reasons. This guarantees that Y 1 > Y 2 • The 
algorithm on the following page implements this theorem. 



5. GENERALIZATIONS 

In all of the processes discussed here except the 
exponential GNEAR(l) the correlations are non-negative and 
geometrically decreasing. A particular scheme for obtaining 
negative correlation is given for the exponential case. 
Another scheme for obtaining alternating correlations which 
are possibly negative and which is broadly applicable is 
given in Gaver and Lewis [1] and in Lawrance and Lewis [5]. 
Another problem is that different types of dependence and 
higher-order Markovian dependence might be encountered in 
data. Schemes for obtaining mixed autoregressive moving 
average exponential sequences where the autoregression has 
order p and the moving average has order q are given in 
Lawrance and Lewis [4]. The mixed exponential process 
TMEAR(l) is easily extended to give a process with this type 
of extended correlation structure. This will be discussed 
elsewhere. 
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ALGORITHM TMEAR1 



This algorithm generates a sample of size N from the TMEAR1 
process. It is based on the generation of arrays of uniform and unit 
exponential numbers. Subroutines UNIFORM and EXPON are assumed to exist 
to generate such arrays. For each TMEAR1 number generated two raw uniform 
random numbers and one exponential number are needed. 



INPUT N, PI1, PI2 , MU1, MU2, ALPHA 



MU PI1 * MU2 + PI2 * MU 2 

B MU1 + MU2 - ALPHA * MU 

A (1 - ALPHA) * MU1 * MU2 
T •«- SQRT(B * B - 4 * A) 

G2 .5 * (B - T) 

G1 *■ .5 * (B + T) 

GO -*• PI2 * MU1 + PI1 * MU2 
El «- (G1 - GO) /G1 - G2) 

R1 = 1/G1 
R2 = 1/G2 

CALL UNIF0RM(U,2*N+1) 

CALL EXPON (E,N+1) 

IF U (2N+1) £ PI1 THEN 

ELSE 



/* Initialization */ 



/* Initialization */ 
/* Generate 2N+1 uniforms */ 

/* Generate N+l unit exponentials */ 
X(0) «- MU1 * E (N+l) 

X(0) -s- MU2 * E (N+l) 



END IF 

DO I 1 to N 

IF U(I) <_ El THEN Z <- R1 *E(I) 

ELSE Z +■ R2 *E (I) 



END IF 

IF U(I+I) £ ALPHA THEN Y ^ (X(I-l) 

ELSE Y 0 



END IF 

X(I) ■«- Y + Z 



END DO 
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